Harnessing Brazilian biodiversity database: identification of flavonoids as potential inhibitors of SARS-CoV-2 main protease using computational approaches and all-atom molecular dynamics simulation

SARS-CoV-2 (Severe Acute Respiratory Syndrome Coronavirus 2) is the etiological agent responsible for the global outbreak of COVID-19 (Coronavirus Disease 2019). The main protease of SARS-CoV-2, Mpro, is a key enzyme that plays a vital role in mediating viral replication and transcription. In this study, a comprehensive computational approach was employed to investigate the binding affinity, selectivity, and stability of natural product candidates as potential new antivirals acting on the viral polyprotein processing mediated by SARS-CoV-2 Mpro. A library of 288 flavonoids extracted from Brazilian biodiversity was screened to select potential Mpro inhibitors. An initial filter based on Lipinski’s rule of five was applied, and 204 compounds that did not violate any of the Lipinski rules were selected. The compounds were then docked into the active site of Mpro using the GOLD program, and the poses were subsequently re-scored using MM-GBSA (Molecular Mechanics Generalized Born Surface Area) binding free energy calculations performed by AmberTools23. The top five flavonoids with the best MM-GBSA binding free energy values were selected for analysis of their interactions with the active site residues of the protein. Next, we conducted a toxicity and drug-likeness analysis, and non-toxic compounds were subjected to molecular dynamics simulation and free energy calculation using the MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) method. It was observed that the five selected flavonoids had lower MM-GBSA binding free energy with Mpro than the co-crystal ligand. Furthermore, these compounds also formed hydrogen bonds with two important residues, Cys145 and Glu166, in the active site of Mpro. Two compounds that passed the drug-likeness filter showed stable conformations during the molecular dynamics simulations. Among these, NuBBE_867 exhibited the best MM-PBSA binding free energy value compared to the crystallographic inhibitor. Therefore, this study suggests that NuBBE_867 could be a potential inhibitor against the main protease of SARS-CoV-2 and may be further examined to confirm our results.


Introduction
In late 2019, a previously unknown zoonotic disease called Coronavirus Disease 2019 (COVID-19) surfaced as a viral respiratory infection in Wuhan, China.Coronavirus 2 (SARS-CoV-2), rapidly spreading worldwide.In March 2020, the World Health Organization (WHO) declared this infectious disease a pandemic (Chauhan, 2020).As of 25 October 2023, there have been 771,549,718 confirmed cases of COVID-19, resulting in 6,974,473 deaths worldwide, according to the WHO (https://covid19.who.int, accessed on 27 October 2023).Although COVID-19 vaccines have been developed remarkably quickly, it remains uncertain how well they will perform in the long term and whether they will be effective against potential future variants of SARS-CoV-2 and other coronaviruses (Heaton, 2020;Jackson et al., 2020;Yu et al., 2020;Panda et al., 2023).With commercially available antiviral treatments proving to be extremely ineffective (Liu et al., 2022;Purohit et al., 2023), the need for the search for potentially more effective antivirals arises.
The primary goal of antiviral drug research for SARS-CoV-2 is to target the viral entry mechanism or viral genome replication (Asselah et al., 2021).SARS-CoV-2 possesses an RNA genome encoding two polyproteins, pp1a and pp1ab, four structural proteins, and various accessory proteins (Qiao et al., 2021).The main protease of SARS-CoV-2 (Mpro, also known as 3CL protease) cleaves the polyproteins at 11 different sites to release mature non-structural proteins (nsps), which play crucial roles in viral RNA replication and immune evasion (Jin et al., 2020;Owen et al., 2021;Wu et al., 2022).Inhibiting Mpro would essentially disrupt viral proliferation throughout the human body (Zhang et al., 2020;Singh et al., 2022).Considering the essential role of Mpro in maturing nsps and SARS-CoV-2 replication, combined with the absence of known human proteases with similar cleavage specificity, makes it an attractive target for the development of antivirals without adverse effects (Qiao et al., 2021;Hu et al., 2022;Liu et al., 2022;Huang et al., 2023).Several SARS-CoV-2 Mpro inhibitors have undergone clinical trials, and nirmatrelvir and ensitrelvir have received regulatory approval (Owen et al., 2021;Mukae et al., 2022;Mukae et al., 2023).However, mutations conferring resistance to nirmatrelvir and ensitrelvir have already been identified in SARS-CoV-2 Mpro (Heilmann et al., 2022), highlighting the importance of seeking new Mpro inhibitors.
Due to the urgent need to find an effective therapy to combat SARS-CoV-2, several studies have emphasized the importance of Natural Products (NPs) as potential antiviral agents (Antonio et al., 2020;Islam et al., 2020;Mani et al., 2020;Wang and Yang, 2020;Rubio-Martínez et al., 2021).Furthermore, some research has focused on the use of NPs as inhibitors of the main protease of SARS-CoV-2 (Mpro) (Aanouz et al., 2021;Ray et al., 2022;Ahamed and Arif, 2023;Purohit et al., 2023).
Among readily accessible Natural Products (NPs), flavonoids stand out, which are naturally occurring polyphenolic compounds found in fruits and vegetables, and have demonstrated significant antiviral activities (Zakaryan et al., 2017;Cataneo et al., 2019;Wang and Yang, 2020).These compounds play a relevant role in enhancing the host's defense system against viral infections, reducing infection, and inhibiting cytokine generation (Jo et al., 2020;Solnier and Fladerer, 2021;Alzaabi et al., 2022).The potential of flavonoids as antiviral agents against SARS-CoV-2 is reported in the literature (Cherrak et al., 2020;Wang and Yang, 2020;Gogoi et al., 2021;Liskova et al., 2021;Owis et al., 2021;Alzaabi et al., 2022;Liu et al., 2022;Chatterjee et al., 2023).Jo and co-workers (Jo et al., 2020) identified a series of flavonoids as effective inhibitors of SARS-CoV-2 replication in in vitro experiments.Furthermore, Mohapatra and colleagues discussed the specific importance of the Cys145 residue in likely covalent interactions with natural products (Mohapatra et al., 2023).Recently, in a SARS-CoV-2 replication assay and molecular dynamics studies, flavonoids were identified as inhibitors of Mpro (Lin et al., 2023).Recent reviews highlight flavonoids with confirmed in vitro anti-SARS-CoV-2 Mpro activity (Solnier and Fladerer, 2021;Antonopoulou et al., 2022).Therefore, exploring flavonoids as potential inhibitors of SARS-CoV-2 Mpro emerges as a promising strategy in the search for alternative therapies based on natural products.
The present study aimed to employ computational approaches to screen flavonoids derived from Brazilian biodiversity with the purpose of identifying potential inhibitors of the main protease (Mpro) of SARS-CoV-2.The results of this study will guide future efforts in identifying and developing Mpro inhibitors based on natural flavonoids.Additional experimental studies will be necessary to validate computational predictions and assess the activity of these compounds in cellular models.
(NuBBEDB) was used in the virtual screening (Valli et al., 2013;Pilon et al., 2017).Subsequently, all compound structures were filtered based on the Lipinski "Rule of 5" (RO5) physicochemical filter (Lipinski et al., 2001) to obtain drug-like molecules.Compounds that did not violate Lipinski's rule were considered for the virtual screening study.
The 3D crystal structure of SARS-CoV-2 Mpro with PDB ID 6M2N complexed with the crystallographic inhibitor Baicalein (3WL) was used for virtual screening (Su et al., 2020).The Chimera software was used for protein preparation, removing the crystallographic ligand 3WL and water molecules.Hydrogen atoms were added, and protonation states at pH 7.0 were assigned to the protein using the H++ server (Gordon et al., 2005).
The GOLD program (Verdonk et al., 2003) was used to perform molecular docking simulations.First, the crystallographic ligand 3WL was re-docked, and the RMSD between the crystallographic ligand and the generated pose was calculated to validate the molecular docking protocol.The docking spatial search docking sphere was defined using the coordinates of the 3WL ligand in the protein's crystallographic structure with a radius of 10 Å using the GoldScore scoring function (Verdonk et al., 2003).After the validation step, natural products from the flavonoid classes were subjected to molecular docking following the same protocol as validation.Subsequently, the poses were rescored according to the workflow proposed by Harutyun Sahakyan, 2021;Sahakyan, 2021), which uses the iPBSA code to reclassify the poses generated in molecular docking through the calculation of binding free energy using the MM-PB (GB) SA methods.
For the calculation of MM-GBSA binding free energy, protein and ligand parameterization was performed first, then the complexes were minimized, and finally, MM-GBSA calculations were performed using AmberTools23 (Case et al., 2023).Ligand charges were calculated using the AM1-BCC method ( (Jakalian et al., 2002).The ff14SB force field was used to describe protein parameters, and the General Amber Force Field (GAFF) for ligands (Wang et al., 2004;Maier et al., 2015).Required input files (coordinates and topologies) with mbondi3 radii were prepared using tLEaP.Minimization was carried out in generalized implicit solvent models for each protein-ligand complex using the sander mechanism (Salomon-Ferrer et al., 2013).For each complex, the maximum number of minimization cycles was set at 1,000, with 500 steps using steepest descent and another 500 steps using the conjugate gradient algorithm.The final snapshot of the minimized protein-ligand complex was used to rescore the poses through MM-GBSA binding free energy calculation.

Molecular dynamics simulation
The compounds that passed the toxicity analysis were then subjected to a molecular dynamics (MD) simulation study to analyze the stability and potential conformational changes of the target protein complexed with flavonoids using Amber22 (Costa et al., 2023).
The structures of the protein and ligands were treated using the Amber ff14SB and the general Amber force field (GAFF), respectively.For MD, the systems were solvated in the tLeap module in a truncated octahedral box with a 12Å edge and filled with TIP3P water (Price and Brooks, 2004).Counter ions were added as needed to neutralize the systems.Subsequently, the systems were minimized in four steps, gradually removing position restraints on atom groups.The minimization process starts with the steepest descent algorithm and switches to the conjugate gradient.Next, slow heating for 500 ps to 300 K at constant volume with position restraints on the solute and unrestricted equilibration for 500 ps at constant pressure were performed.Using a collision frequency of 2 cm −1 and linked to a Langevin thermostat, the temperature was maintained at 300 K.The SHAKE algorithm (Ryckaert et al., 1977) and Particle Mesh Ewald (PME) (Darden et al., 1993) were used to limit bond lengths involving hydrogen atoms, and a 10 Å cutoff was applied for non-bonded interactions.Finally, a 100 ns production run was conducted without positional restraints at a constant temperature of 300 K.

Binding free energy calculations
To describe the binding affinity of the natural products and the crystallographic inhibitor with the Mpro protein, we performed binding free energy calculations using the MM/PBSA method (Genheden and Ryde, 2015) available in the AmberTools23 package (Case et al., 2023), while the mathematical description was published in previous studies (Costa et al., 2022).The last 10 ns of the MD simulation trajectories were used to calculate the binding free energy and the decomposition of the binding free energy.

Toxicity analysis
The five flavonoids were analyzed for different types of toxicities, including tumorigenic, mutagenic, reproductive effectiveness, irritant, and drug-likeness using ORISIS Data Warrior 5.2.1 (Sander et al., 2015).

Virtual screening and molecular docking
Although some natural products (NPs) have bioactivity or druglike activity, many NPs have low solubility or chemical instability, making their transformation into parenteral drugs challenging.Additionally, the complex structures of many NPs result in high molecular weights, which are likely to have negative effects on intestinal absorption (Li et al., 2019).
In this study, the bioavailability of the 288 flavonoids was evaluated based on their physicochemical properties using the Lipinski's rule.This step was carried out within the NuBBEDB.It was found that 204 flavonoids did not violate any of the Lipinski's rules and were therefore selected for the virtual screening phase.
Before the virtual screening of the flavonoids, the molecular docking protocol was validated by redocking the co-crystallized ligand 3WL into the active site of Mpro.It was found that the redocked conformation of the ligand perfectly overlapped with the co-crystallized ligand, with an RMSD value of 0.636 Å (Supplementary Figure S2).This result suggests that the GOLD program exhibits satisfactory accuracy in repositioning the 3WL ligand within the active site of Mpro.The docking protocol is considered satisfactory when the RMSD value between the docking pose and the crystallographic ligand pose is less than 2.0 Å (Wang et al., 2016;Costa et al., 2023).Figure 1 illustrates the overlapping alignment of the redocked pose and the crystallographic ligand.
The validated docking protocol was subsequently used for virtual screening based on molecular docking of 204 flavonoids extracted from the Brazilian biodiversity available in NuBBEDB (Pilon et al., 2017).After being scored by the GoldScore scoring function in the GOLD program, the poses were rescored using the method of binding free energy (ΔGbind) calculation MM/GBSA using AmberTools23.
Virtual screening has become a widely applied method in earlystage enzymatic inhibitor discovery projects.The success in predicting hit candidates is closely related to the scoring function applied.However, docking scores have limitations as they often show a weak correlation with experimental results, with only a small fraction of hits exhibiting activity in in vitro validation assays (Rastelli and Pinzi, 2019).Because of this, in our study, we applied the rescoring method using MM/GBSA binding free energy calculations as it is a more accurate method for the selection of the best molecules in virtual screening studies (Rastelli and Pinzi, 2019;Poli et al., 2020;Sahakyan, 2021).
After molecular docking and rescoring of the poses, the top five ranked flavonoids based on the MM/GBSA values were selected, and the interactions of these structures with the amino acid residues in the active site of the Mpro protein were analyzed and compared with those of the crystallographic ligand 3WL, used as a reference compound in this study.
Table 1 shows  S1 shows the ΔGbind values obtained by the MM/GBSA method for all flavonoids used in this study.
The analysis of the new MMGBSA scores, used to avoid false predictions, revealed that the five selected compounds have better binding affinities with the Mpro enzyme, ranging from −51.57 to −47.08 kcal mol −1 , when compared to the reference inhibitor (3WL: −31.26 kcal mol −1 ), as shown in Table 1.
We analyzed the binding modes of these compounds to the active site of the Mpro enzyme.The interactions made by the compounds and the inhibitor 3WL are shown in Figure 2. Our results demonstrated that the selected compounds formed an equal or greater number of hydrogen bonds with the target protein compared to the crystallographic inhibitor 3WL.This explains the higher binding affinity of these flavonoids to the Mpro protein than 3WL (Gogoi et al., 2021).
The analysis of the residues involved in the interactions between the compounds and Mpro showed that NuBBE_ 867 interacted with the residues Ser144, Cys145, Glu166, Gln192, Gln149, and Thr190.NuBBE_1310 formed interactions with the residues Cys145 and Glu166.NuBBE_ 1884 interacted with Gly143, Ser144, Cys145, and Glu166.NuBBE_1890 interacted with the residues Gly143, Cys145, and Glu166, and NuBBE_2328 interacted with Glu166 and Thr190.All of these residues are part of the protein's active site and are responsible for stabilizing these ligands (Su et al., 2020).
The flavonoids NuBBE_867, NuBBE_1310, NuBBE_1884, and NuBBE_1890 form hydrogen bonds with the important catalytic residue Cys145, which is part of the S1 binding pocket present in the catalytic domains I and II of the SARS-CoV-2 Mpro.Along with residue His41, they form the catalytic dyad within the substrate-binding site that actively participates in the catalytic function of the protein.Therefore, binding of these flavonoids to this residue may reduce the catalytic activities of Mpro, ultimately leading to a reduction in viral replication (Gogoi et al., 2021).
Another important residue that interacts with all selected flavonoids is Glu166, a key amino acid involved in the dimerization of the Mpro protein and assists in creating the S1 substrate binding pocket (Goyal and Goyal, 2020;Silvestrini et al., 2021).The formation of a hydrogen bond with the Glu166 residue is reported in the study by Duong and Nguyen (Duong and Nguyen, 2023) as one of the key interactions of noncovalent Mpro inhibitors.It is interesting to note that the crystallographic inhibitor 3WL, used as a reference in this study, also formed a hydrogen bond with the Glu166 residue, which corroborates our results (Su et al., 2020).
The residues Cys145 and Glu166 have been recently reported in molecular docking studies of inhibitors involving natural flavonoids as key residues in the binding process at the Mpro site (Cherrak et al., 2020;Su et al., 2020;Lin et al., 2023), and also with different compounds, corroborating with our results (Gahlawat et al., 2020;Kumar et al., 2020;Gogoi et al., 2021;Ngo et al., 2021;Patel et al., 2021;Antonopoulou et al., 2022;Salarizadeh et al., 2022;Ahamed and Arif, 2023;Rampogu et al., 2023).It is important to mention that all ligands fit perfectly into the protein's binding pocket, revealing that all ligands form a stable complex with the target receptor.Validation of molecular docking protocols using the GOLD program for the crystal structure of the COVID-19 main protease (Mpro) in complex with 3WL inhibitor.Blue is the co-crystal ligand and green is the docking pose.
Frontiers in Chemistry frontiersin.org04 da Rocha et al.
10.3389/fchem.2024.1336001TABLE 1 Binding free energies of the top 5 rescored hits along with the co-crystallized ligand 3WL obtained using the MM/GBSA method.
The compounds with best scores Binding free energy (kcal mol

Toxicity analysis
The five flavonoids were subjected to drug-likeness evaluation and assessment of mutagenic, tumorigenic, reproductive effective, and irritant parameters using ORISIS Data Warrior.It was observed that two flavonoids, NuBBE_867 and NuBBE_1890, showed no toxicity based on the toxicity parameters used in the study.On the other hand, NuBBE_1310 exhibited mutagenic and tumorigenic properties, while NuBBE_2328 displayed mutagenic, tumorigenic, and irritant properties.NuBBE_1884 showed a negative drug-likeness probability.The results of toxicity prediction and druglikeness property analysis are presented in Table 2. Based on the data obtained from ORISIS Data Warrior (Sander et al., 2015), compounds with higher or positive drug-likeness probability values are considered good drug candidates.Since compounds NuBBE_1310 and NuBBE_2328 exhibited toxic effects, and NuBBE_1884 had a negative drug-likeness probability, they were not considered for further analysis.The non-toxic compounds, NuBBE_867 and NuBBE_1890, underwent molecular dynamics simulation studies.
TABLE 1 (Continued) Binding free energies of the top 5 rescored hits along with the co-crystallized ligand 3WL obtained using the MM/GBSA method.

The compounds with best scores
Binding free energy (kcal mol

Molecular dynamics simulation
Molecular Dynamics (MD) simulation was performed to assess the flexibility, and stability of the protein-ligand complex, and the interactions of ligands with the SARS-CoV-2 Mpro protein.In this study, MD simulations were conducted for the Mpro-NuBBE_867 and Mpro-NuBBE_1890 complexes, selected based on toxicity and drug-likeness analyses, as well as the Mpro-3WL (reference inhibitor) complex and the unbound form of the protein (Apo).To predict the dynamic behavior and structural changes in the  complexes, MD simulations were carried out over a 100 ns and the resulting trajectories were used to investigate complex RMSD, RMSF, and the number of hydrogen bonds.
RMSD analysis was performed to monitor the conformational stability of the ligand-protein complexes (Mpro).The RMSD plots for the ligand-protein complexes (Mpro-NuBBE_867, Mpro_ NuBBE_1890, and Mpro-3WL) and the Apo form are shown in Figure 3.The stability of the complexes is directly related to the RMSD values, with lower deviation indicating greater structural stability.The range that justifies complex stability is approximately ~3 Å (Rafique et al., 2023).Our results demonstrate that the average RMSD values for all complexes fall within the range of 1.47-2.14Å, indicating that all complexes remain stable during the 100 ns simulation trajectory (Figure 3).
The average RMSD values for the Apo protein were 1.89 Å, 2.04 Å for the reference complex (Mpro-3WL), 1.47 Å for Mpro-NuBBE_867, and 2.14 Å for Mpro-NuBBE_1890.The RMSD values revealed that the natural products (NPs) NuBBE_867 enhanced the structural stability of the Mpro protein (Figure 3), inhibiting its active site when compared to the crystallographic inhibitor used as a reference.On the other hand, NuBBE_1890 showed an average RMSD value higher than the Apo protein and the reference inhibitor, indicating that the presence of the ligand influenced the stability of the Mpro enzyme and altered its dynamic behavior.The results of RMSD values showed that the inhibitors formed stable bonds with the active site residues of the protein.
One essential component in determining the stability of a protein-ligand complex is the Root Mean Square Fluctuation RMSD of the Apo form of Mpro and the Mpro complexes with 3WL, NuBBE_867, and NuBBE_1890 over a 100 ns MD simulation time.
Residues play a crucial role in forming a stable and strong binding complex between a protein and a ligand (Surti et al., 2020;Singh et al., 2022).The RMSF for the Apo protein and all simulated complexes are displayed in Figure 4.The RMSF of 3WL showed slightly higher fluctuations when compared to the other simulated systems, suggesting that this system is less stable than the others, a hypothesis that is later supported by the free energy calculations.In general, regions with lower fluctuations included the active site regions of Mpro, including Glu166, Cys145, Asn142, Phe140, Met165, Met49, His41, and Leu141, which comprise more stable regions.In regions where the simulated ligand interacts with amino acid residues of the protein, lower fluctuation values indicate greater stability in those interactions.Particularly noteworthy are the regions of the Ser46, Met49, Cys145, and Gln189 residues.
The fluctuations of the binding site residues Glu166, Cys145, Asn142, Phe140, Met165, Met49, His41, and Leu141 were lower in all Mpro complexes.Small fluctuations in these binding site residues indicated high stability of the Mpro complexes.Furthermore, the data suggested that the residues involved in the binding were responsible for the stable RMSD of the Mpro complexes.

Hydrogen bonding
The analysis of hydrogen bonds is an essential parameter that provides information about the binding affinity of drug candidates to a protein.Therefore, the presence of a significant number of H-bonds indicates a strong interaction between a ligand-protein complex.It is also responsible for drug metabolism and specificity (Rafique et al., 2023).
We investigated the hydrogen interactions of the flavonoids NuBBE_867 and NuBBE_1890 with the Mpro protein's site residues and compared them to the crystallographic inhibitor 3WL during a 100 ns MD simulation (Figure 5).
Our results show that NuBBE_867 formed five hydrogen bonds with the active site of Mpro throughout the simulation, a higher number of interactions compared to NuBBE_1890 and the crystallographic inhibitor 3WL, which had three and two interactions throughout the simulation, respectively.
During the simulation period, NuBBE_867 interacted with the site residues Asp187, Glu166, Ser144, Gly143, and Cys145 with occupancies of 54.19%, 43.06%, 12.28%, 12.24%, and 11.48%, respectively.NuBBE_1890 interacted with the residues Asn142, Thr24, and Glu166 with occupancies of 17.82%, 7.89%, and 6.29%, respectively.In contrast, the crystallographic inhibitor interacted with the residues Arg188 and Glu166 with occupancies of 21.17% and 11.10%, respectively.The formation of multiple hydrogen bonds (H-bonds) between the NuBBE_867 inhibitor and the Mpro protein contributed to the greater stabilization of the structural complex, as shown in the RMSD analysis.The residues Gly143, Ser144, and Cys145 that interacted with NuBBE_ 867 constitute the "oxyanion hole" of this cysteine protease and play a fundamental role in catalyzing the hydrolysis of protein substrates, in line with studies on Mpro inhibitors by flavonoids (Lin et al., 2023).It is important to note that the interactions with the residues Cys145 and Glu166 are consistent with molecular docking.The Cys25 residue is essential for the catalytic function of Mpro, while the Glu166 residue is crucial for Mpro dimerization and the formation of the substrate binding pocket (Cherrak et al., 2020).In a previous study, Yoshino and colleagues ((Yoshino et al., 2020) conducted molecular dynamics simulations involving different inhibitors complexed with Mpro to verify the necessary interactions with the active site residues of Mpro for considering a compound as an inhibitor.They revealed that residues Glu166, Gly143, Ser144, and Cys145 were among the main interacting residues, which corroborates our results.Similar interactions have also been observed in previous studies involving natural flavonoids, further reinforcing our findings (Cherrak et al., 2020;Su et al., 2020;Lin et al., 2023).Our results demonstrate that NuBBE_867 effectively interacted with the binding pocket of the Mpro protein, is highly stable, and is expected to be a potential candidate as an inhibitor of the main protease of SARS-CoV-2.

MM/PBSA binding free energy
To confirm the inhibitory capacity of the ligands for the Mpro protein in the ligand-protein complexes, binding free energy calculations were performed.The calculated results of MM-PBSA binding free energy are presented in Table 3.
The MM-PBSA calculations results for NuBBE_867, NuBBE_ 1890, and the inhibitor 3WL are shown in Table 3. From the results, the ΔGbind of the Mpro-3WL complex was −15.97 kcal/mol.The ΔGbind of Mpro-NuBBE_867 and Mpro-NuBBE_1890 was −30.04 kcal/mol and −24.92 kcal/mol, respectively.It was observed from the MM-PBSA analysis that the complexes formed between the flavonoids and Mpro had a lower ΔGbind than the reference complex, Mpro-3WL.This indicates the formation of stable complexes with higher binding affinity for these flavonoids in the active site of Mpro.It's worth noting that NuBBE_ 867 exhibited the highest binding affinity to Mpro, suggesting a more potent inhibitor of Mpro when compared to the 3WL inhibitor.
Notably, these results are in line with the findings in the analysis of hydrogen bonds, which show that NuBBE_867 formed interactions with key residues in the active site of Mpro, explaining the higher affinity with Mpro.The MM-PBSA method has been effectively used in previous computational studies to estimate and select more potent inhibitors of the Mpro protein from natural sources, further supporting our results (Ahmad et al., 2021;Gogoi et al., 2021;Sharma et al., 2022;Lin et al., 2023;Rafique et al., 2023).
Additionally, energy terms such as van der Waals energy, electrostatic energy, SASA energy, and polar solvation energy are crucial contributors to the total binding free energy of the complexes.We can see in Table 3 that van der Waals, electrostatic, and SASA energies contributed significantly to the total binding energies of the complexes.On the other hand, polar solvation energy had an unfavorable impact on the binding energy.
To delineate the specific contribution of each amino acid residue to ΔGbind, we conducted an energy decomposition analysis for each residue within the protein.Residues contributing to binding with values of −1.00 kcal/mol or lower were deemed significant for the binding process.The energetic contributions of each residue in the Mpro-NuBBE_867 and Mpro-NuBBE_1890 complexes were compared to those in the complexes formed by the reference compounds Mpro-3WL, as illustrated in Figure 6.
The energy decomposition analysis by residues unveiled that NuBBE_867 engages with His41 (ΔGtotal = −1.77)and Cys145 (ΔGtotal = −1.83),both of which constitute the catalytic dyad within this cysteine protease.Additionally, it forms interactions with Gly143 (ΔGtotal = −1.83), a residue that, alongside Cys145, plays a pivotal role in the oxyanion hole responsible for ligand stabilization.
Figure 6 depicts a greater number of residues with substantial contributions to the flavonoid NuBBE_867 in comparison to the reference ligands 3WL and NuBBE_1890.This observation elucidates the heightened binding affinity of this ligand with the Mpro enzyme, as indicated in Table 3, thus solidifying the validity of our findings.
Our in silico results propose that the molecule NuBBE_867 holds significant promise as an inhibitor of the Mpro enzyme.

Conclusion
We have shown that Flavonoids may be effective to inhibit the SARS-CoV-2 Mpro protein.Molecular docking results indicate that all the compounds exhibit reasonably strong binding energies, ranging from −5.9 to −7.3 kcal/mol, and lower predicted inhibition constants compared to two FDA-approved antiviral drugs, remdesivir and nirmatrelvir, currently used in SARS-CoV-   Graphical representation of the interaction energy per residue.

FIGURE 4
FIGURE 4Root Mean Square Fluctuations of the Apo protein and Mpro-ligand complexes' residues.

FIGURE 5
FIGURE 5Number of hydrogen bonds between the ligands and Mpro protein during the 100 ns MD simulation time.
2 treatment.Subsequent MD simulations confirm the stability of the ligand-protein complexes.Thus, RMSD, RMSF, radius of gyration, and hydrogen bond interactions demonstrate the maintenance of structural stability and productive interactions.In conclusion, based on the collective findings of this study, it is evident that these compounds hold promise as candidates for the development of effective medications to inhibit the SARS-CoV-2 Mpro protein.Further in vitro and in vivo experiments are warranted to validate and advance these results, potentially contributing to the development of novel antiviral agents.

TABLE 2
Toxicity and drug-likeness analysis.

TABLE 3
MM-PBSA binding free energies (ΔG bind ) and their components for the complexes under study.All values are reported in kcal.mol−1 .